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. Abstract We use full available array of radial velocity data, including recently 

published HARPS and Keck observatory sets, to characterize the orbital configu- 
ration of the planetary system orbiting GJ876. First, we propose and describe in 
detail a fast method to fit perturbed orbital configuration, based on the integra- 
tion of the sensitivity equations inferred by the equations of the original iV-body 
problem. Further, we find that it is unsatisfactory to treat the available radial 
p ^ velocity data for GJ876 in the traditional white noise model, because the actual 

noise appears autocorrelated (and demonstrates non- white frequency spectrum). 
The time scale of this correlation is about a few days, and the contribution of 
the correlated noise is about 2 m/s (i.e., similar to the level of internal errors in 
the Keck data) . We propose a variation of the maximum- likelihood algorithm to 
estimate the orbital configuration of the system, taking into account the red noise 
effects. We show, in particular, that the non-zero orbital eccentricity of the in- 
nermost planet d, obtained in previous studies, is likely a result of misinterpreted 
red noise in the data. In addition to offsets in some orbital parameters, the red 
noise also makes the fit uncertainties systematically underestimated (while they 
are treated in the traditional white noise model). Also, we show that the orbital 
eccentricity of the outermost planet is actually ill-determined, although bounded 
by ~ 0.2. Finally, we investigate possible orbital non-coplanarity of the system, 
and limit the mutual inclination between the planets b and c orbits by 5° — 15°, 
| depending on the angular position of the mutual orbital nodes. 

■ Keywords extrasolar planets ■ radial velocity • red noise ■ mean-motion resonance 

o 
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1 Introduction 

The first extrasolar pla net in the system or bitin g the red dwarf GJ 876 was de- 
tected independently bv lDelfosse et all ( 199ct) and lMarcv et al.l ( 1998h on the basis 
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of high-precision radial velocity (hereafter RV) Doppler measurements, acquired 
at ELODIE+CORALIE and Keck HIRES s pectrographs. Soon after this, the sec- 
ond planetary companion was discovered bv lMarcv et al. (2001). That discoveries 
showed that the exoplanetary system of GJ876 is an extraordinary object. The two 
planets b and c were massive giants (having minimum masses m sin i of roughly 2 
and 0.6 of Jupiter masses) orbiting in short period orbits (approximately 60 and 
30 days). Therefore, it became clear that the planets orbital dynamics should be 
significantly affected by the 2/1 mean-motion resonance (hereafter MMR). The 
resonance itself is not yet very astonishing, because now we already know many 
extrasolar planets in various MMRs. The extraordinarity of the GJ876 system 
comes from the fact that interplanetary gravitational perturbations were directly 
detected in the observed radial velocity data, i . e. they reveal themselves on th e 
observational time scale ( Lauehlin fc Chambers!. [200 luRi vera &: Lissaueit l200lh . 

Still, no other extrasolar planetary system around a main sequence star is 
known to demonstrate so clear and well-measurable signatures of planetary per- 
turbations in its radial velocity time series (although, there are some exoplanetary 
systems where such perturbations are suspected to be non-negligible). A few fa- 
voring factors meet each other in GJ876: large masses of the planets b and c, small 
star mass (hence especially large planet/star mass ratios), rather short orbital pe- 
riods (thus shorter perturbations time scale), and, of course, the low-order MMR. 
The dynamical perturbations reveal themselves mainly in the form of secular cir- 
culation of the planetary apsidal lines, triggering slow change in the non-sinusoidal 
shape of the RV oscillations being observed. The notion "secular" , however, looks 
somewhat odd here, since the period of this "secular" circulation is only about 
10 yrs. During the whole observation term since 1998 till today, these apsidal lines 
completed roughly a single revolution each. 

This dynamical effect complicates the procedure of the RV data analysis, but 
in exchange it enables us to determine (solely from the RV data) the inclination 
of the system to the sky tangent plane. This allows us to determine, instead of 
the minimum planetary masses msini, the true masses m, which normally remain 
unconstra ined in the RV exop lanet detections. 

Later, iRivera et al.l (|2005| ) reported the discovery of the third planet in the 
system, based on further Keck RV observations of GJ876. This planet (d) pos- 
sesses very low mass of ~ 7.5 Earth masses and a very sh ort orbital period of 
approximately 2 days. In the paper ( Bean fc Seifahrtl 120091) . the question of or- 
bital non-coplanarity of this system was studied extensively, and these authors 
gave an estimation of ~ 5° for the mutual orbital inclinatio n of the plan ets b 
and c. This result was based on th e Keck RV data from ( Rivera et all 20051) and 
on the HST astrometry data from ( Benedict et al.Ll2002h Recently. ICorreia et al.l 
(2010) reanalyzed these old Keck data adding to them a new array of very accurate 
RV measurements, obtained by the HARPS spectrograph (ESO). He presented im- 
proved three-planet orbital fits, which generally agree with the fits by Rivera et al. 
([20051) . 

Finally, in the very recent paper ( Rivera et al.l[201ol) the discovery of the fourth 
Uranus-mass planet e was reported. This last planet is also remarkable, because it 
forms a Laplace three-planet resonance with two other giant planets, so that the 
ratio of the periods P c : P b : P e is close to 1 : 2 : 4. 

It is interesting, however, that the RV si gnature of this fourth planet had been 
found in the old Keck data already, e.g. in ( Baluevl [20 08b) . It is notable that the 
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fourth pl anet paramet e rs fro m ( Bafuevl . l2008bJ) are surprisingly close to the ones 



given bv iRivera et al.l ( 2010l) . including even the mean longitude and the small 



orbital eccentricity. This o rbital configuration appears pretty stable at long time 
scales. IRivera et al.l ( 2010l) also mention that they were seeing the signs of the 



fourth planet in their data since 2005, but until 2010 they could not find a stable 
four-planet solution, due to the large plan et e eccentric i ty. Th e reason of such 
difference probably comes from the fact that IRivera et al. (2010) did not take into 



ac count the annua l errors of a few m/s, which were detected in the old Keck data 
bv lBaluevi (|2008bh . Such annual errors may appear rel atively frequently in the RV 



planet search surveys, as discussed in ( Baluevl . 12009) . They may have different 



physical sources; in case of GJ876 they were probably caused by some err ors in 
the data reduction pipeline, since new Keck data from IRivera et all ( 2010|) look 



free from such errors. In fact, a careful error analysis of the RV data could enable 
the robust detection of the fourth planet several years ago already. 



Since the quality of the Keck data is considerably improved in (Rivera et al 



2010) , and the completely new very accurate RV data is now published (jCorreia et al.l 
20101 ). we conclude that it is time to perform such detailed analysis now. Accurate 
and unbiased values of planetary masses and orbital parameters in so remarkable 
system as GJ876 represent a huge interest for the celestial mechanics studies, as 
well as for more astrophysical branches like the tidal star-planet interaction stud- 
ies (note the "hot super-earth" GJ876 d) or for the planet formation theories. In 
particular, it is important to have some reliable estimations or limits of the orbital 
non-coplanarity of such system, or to know how large the orbital eccentricity of 
the "hot Earth" GJ876 d actually can be. 

The plan of the paper is as follows. First, in Sect. [2] we describe in detail the fast 
approach of fitting the dynamically perturbed planetary configurations on the basis 
of the star RV time series. We also explain the conventions we adopt to determine 
the reference osculating orbital parameters, and describe the numerical integrator 
we use in this study. In Sect. [3] we describe in more details the RV data we use in 
the paper and also give some preliminary analysis results, obtained during direct 
fitting of these data. In Sect.|H we find that these RV data cannot be modeled in 
a traditional way, because they contain significant fraction of autocorrelated (non- 
white) noise. Also, we propose here a modified RV fitting algorithm, which allows to 
take such correlated RV errors into account properly. In Sect. [5] we study the planet 
e orbital parameters and show that its eccentricity is still ill-determined although 
bounded by 0.2 from the upper side. In Sect. [6] we give final estimations of all 
planetary parameters, taking into account the effects of the RV noise correlateness 
and of the bad determinability of the planet e eccentricity. In Sect. [71 we study the 
(non-)coplanarity of the system in view of currently available RV data. Finally, 
Sect. [8] describes briefly the dynamical evolution of the admissible GJ876 orbital 
configurations. 



2 Fitting graviationally perturbed orbits 

2.1 Evaluating radial velocity function 



Before we proceed to the investigation of the GJ876 planetary system itself, 
we describe the algorithm we used to perform iV-body fitting of RV data. Al- 
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though such iV-body fitting; have been already carried out in several previous 
works (|Laughlin fc Chamber j, l200ll H lvera fc LissauerikoOlh . but they omi t im- 



porta nt details of the algorithms used. We may highlihgt a very recent paper ( Pal, 
20IfJ ). which introduces a method of orbital fitting, which is similar (in spirit) to 
the metho d that we use here. However, many ideas of the methods developed by 



Pal (|2010h look pretty different, so we still need to describe our method in detail 



here. The algorithm that we use represents a variation of the popular non-linear 
Levenberg-Marquardt x 2 minimization algorithm with a minor modificat ion nec- 
essar y to take into account the "RV jitter" phenomenon, as discussed in (jBaluevl 
l2009h and below. This algorithm is an iterative gradient method, which requires 
to evaluate (on different iteration stages) the fittable model (RV curve model, in 
our case) and its partial derivatives over all free parameters. In case of GJ876, 
the two latter subtasks are non-trivial, since we need to take into account the full 
dynamical model of planetary perturbations. 

Let us write down the equations of the iV-body problem in the astrocentric 
coordinate system: 

» 3=1. .AT \ 1 J U 3 J 

Here /ij are the ratios of the planet masses m, to the star mass M*, and k 2 = GM* 
is the analogue of the Gauss gravitational constant. Note that further we may call 
fj,i just "planetary masses" for shortness, since usually it should not introduce any 
misunderstandings. We assume that z axis is directed along the observer's line of 
sight, and two other axes are directed arbitrarily in the tangent sky planeQ 

The initial conditions for the system ((T|) may be expressed via some set of 6jV 
osculating orbital elements forming M vectors -Ki using the well-known functions 
of the Keplerian motion. However, there is no unique definition of the osculating 
Keplerian elements. Their set depends not only on our choice of the parametriza- 
tion, but also on the non-unique splitting of the forces in ([!]) into the unperturbed 
Keplerian and perturbational parts. For instance, we can assume that the Keple- 
rian part of the force is — fc 2 ri/r 3 , and everything else in the right hand side of ([1]) 
is the perturbation. In such a case, the initial planetary coordinates and velocities 
could be expressed via corresponding osculating orbital elements as 

ri| t=0 = fc^Krfa), v;| t=0 = fc5K fa), (2) 

where the functions K r and K. v should represent the solution of the standard equa- 
tion of the Keplerian motion r = — r/r 3 at t = 0. Each vector tt incorporates six 
osculating orbital elements of the planets: orbital period P, eccentricity e, percen- 
ter argument w, mean argument of latitude ip (mean anomaly + w), inclination to 
the sky plane i, longitude of the ascending node Q. The specific of the exoplane- 
tary orbit determination problem sets a few pecularities of these parameters that 
we need to highlight: 



1 Note that the term "iV-body" , which we use in the paper frequently, does not imply that 
the number of the bodies in the system is denoted by N. Instead, J\f stands for the number of 
interacting planets, N for the number of RV measurements, and 'W-body" must be understood 
as a solid term. 
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1. We remind the reade r that the angle uj i s tr aditionally counted from the de- 
scending orbital node ( Ferraz-Mello et all I2005L sect. 2.2). Historically, this de- 



viation from the classical definition of uj occured because exoplanet researchers 
traditionally model the host star's RV wobble using the Keplerian formula 
V = K(cos(uj + v) + ecoscj), which actually expresses the RV of the orbit- 
ing planet itself. This shifts the derived value of uj by n, due to the RV sign 
change. Alternatively, we may think that uj refers to the star's orbit around the 
common barycenter, although such treatment becomes a bit vague for multi- 
planet systems. This pecularity of the angle uj is insignificant for the most of the 
exoplanets, since the positions of the planetary orbital nodes usually remain 
unconstrainable anyway. It may become important for GJ876, however. 

2. Consequently, the angle ip is also counted from the descending node. We use ip 
instead of the full mean longitude A = ft + ip + iv, since the absolute positions of 
the orbital nodes still remain unconstrained until the astrometric observations 
are used. Radial velocities can only constrain the differences between J7j. For 
coplanar configurations, dealing with ip and oj is equivalent to dealing with A 
and the pericenter longitude -uj = fl + uj + ty, since all Qi are equal to the same 
constant value. 

3. In the equations ([2]), the dependence on the actual k (in fact, on M*) is ex- 
tracted separately in the external factors fc 2 / 3 . Such dependence on k appears 
when 7r contains planetary orbital period P (or e.g. mean motion) as a primary 
parameter, and the semi-major axis is treated as only a derived one. We use 
such "period-oriented" choice because it is the planetary orbital period (or the 
period of corresponding RV oscillation) which is drawn from the observations 
directly, rather than the orbital semi-major axis. If we chose the semi-major 
axis as a primary parameter, the factors in (J2j) would be 1 and k. 

The astrocentric coordinate system makes the TV-body equations ([TJ) simple 
for numerical integration, but it is not well suitable to reference orbital elements. 
This is mainly because of relatively significant dependence of the planetary orbital 
per iods on the choice of the coordinate system, which has been already discussed 
in (jLissauer fc Riveral l200lfc iLee fc Pealel. 120031 : iBaluevl l2008cf ). As it was noted 



in these works, the osculating orbital periods (or mean motions) referenced in the 
Jacobi coordinate system generally match better the apparent orbital periods (or 
apparent average drift of the mean longitude). In practice that would mean that, 
e.g., we may freely experiment with turning on/off the gravitational interactions of 
any planet without the need to manually adjust its orbital period when switching 
from one model to another. Otherwise, we should remember that there may be 
a significant offset between the best fitting orbital period in the Keplerian RV 
model framework (which implies apparent orbital period) and the full iV-body 
framework (the osculating period). Although such offsets are small (C(/i,j)), they 
may be comparable to or even significantly exceed the relative statistical error 
of the respective periods (e.g. ~ 10~ 5 in case of the innermost planet of GJ876). 
In such a case, the orbital fitting is slowed down and eventually may even fail 
to converge to the correct period value (it may be attracted by a spurious close 
alias or just noisy neighboring period). The choice of the Jacobi system usually 
eliminates these troubles, if we also assume that each osculating Keplerian orbit 
refers to a fictitious central mass incorporating the star's and this planet's mass 
and masses of all inferior planets. 
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According to (Baluev, 2008c, appendix A), one way to define the unperturbed 



part of the iV-body Hamiltonian in the Jacobi coordinates is 

.EL 

1[L 



ll , 2 i 

rr Ji Pi k 1 I V~* lr>\ 

H Kep ,i = —^--li-i—f-, 7i = l + 2^W- ( 3 ) 



We adopt this (non-unique) way to split the Hamiltonian into the Keplerian and 
perturbational parts, because the corresponding equation of the Keplerian motion 



involves the value of the central mass instead of M*. This is exactly what 

we need. The initial conditions for the Jacobi coordinates and velocities should 
represent a solution of the latter equation and therefore should look like 

- - 

r i| i=Q = k'i 3 K r (n t ), v i| t=0 = ki 3 K v (iri) (5) 

These formulae differ from ([2]) only by the coefficient k, which now involves plan- 
etary masses too. Note that in case of a highly hierarchical system with negligible 
mutual perturbations between planets, such definition of the osculating orbital 
periods would make them infinitesimally close to the apparent revolution periods. 
The exoplanetary systems are not always hierarchical, but in practice (in partic- 
ular, for GJ876) the offsets between apparent and so-defined osculating planet 
periods in -Ki usually remain satisfactory. 

After that, we can transform the Jacobi vectors ([5]) to the astrocentric ones 
using the formulae 

ri = Tj(ri, . ..,Vi,tix,.. .,/J.i-i) = Ti+ } — rj, 

3=1 

v l = Tj(vi, . . . ,v-,^ii, . . (6) 

The formulae ([5]) and ([6]) jointly express the Cartesian initial conditions in the 
astrocentric system via the osculating orbital elements -Ki defined in the Jacobi 
system, and via the planet masses m. For shortness, we can write down these 
compound dependences as 

r i\t=0 = ^U(k,TTl, . . . ,7Ti, fll, . . . ,fli), 

Vi\ t=0 = Vj(fc,7Ti,. ..,7r,,/xi,. (7) 

Note that these functions also involve the parameter k, which depends on the star 
mass A/*, which we assume is a priori known from astrophysical models of its 
spectrum. Following ICorreia et aTl (|2010l) , we assume Af* = O.334M0 throughout 
the paper. 

Having the astrocentric initial conditions from (JT]), we can integrate (fTJ) to 
the desired time. After that, we will have astrocentric planetary positions and 
velocities. We prefer to integrate the astrocentric motion equations (JXJ) , since the 
analogous equations in the Jacobi system would be significantly more complicated 
and thus would slow the algorithm down. After the integration, the barycentric 
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velocity vector of the star can be obviously expressed via the planetary astrocentric 
velocities as 




Finally, all stages explained in this section allow us to determine the model RV 
curve, based on the planetary orbital parameters -k^ and masses fii as fittable 
free variables. The total observable radial velocity of the star incorporates also a 
constant radial velocity (we denote it as co), caused by the motion of the planetary 
system barycenter. Moreover, since the RV measurements available for GJ876 are 
basically relative, we should also introduce different fittable values of co for the 
data series obtained at different instruments. 



2.2 Using variational equations 



Probably the most important stage of the Levenberg-Marquardt algorithm (as well 
as of any other gradient optimization method) is evaluation of the partial deriva- 
tives of the data model over the free parameters. In problems with complicated 
model function (like iV-body RV fitting) most of the computational time is spent 
during evaluation of these derivatives. The simple way to evaluate them is just to 
use numerical differentiation of the original RV model (which should be calculated 
using TV-body integration) . This way is easy to implement algorithmically, but this 
leads to a very bad speed/error ratio, since the error of numerical differentiation is 
considerably larger than the error of the original function. A better way to evaluate 
these derivatives is to integrate variational equations of the TV-body problem (in 
the literature on non-linear optimi zation they a re also called sensitivity equations, 
see e.g. chapter 8 in the book by iBardl \l97^i ). It allows to obtain roughly the 
same accuracy of derivatives as in the original function with the same integration 
step and number of equations to integrate. Planetary positional vectors and 
their velocities Vj = dvi/dt depend on the time, on the planetary masses fn and on 
the orbital parameters -Ki. Differentiation of ([I]) over iTi and fn yields the following 
linear ODE system for partial derivatives of : 



1 <T dr 



1 d 2 drj 
k 2 dt 2 dfij 



Or; 



+ £ 



k^i 



drj 



fc=l..jV 

k^i 



{j %) \dfij dfij) {j) dtij. 



n(a) 



a I — 3a ® a 



l,i = 3 



(9) 



Here 5y is the Kronecker delta, I is the identity matrix, <g> denotes the dyadic prod- 
uct of vectors, and thus l~l(a) is a 3 x 3 matrix having elements (<5jj — "iaiaj/a 2 )/^. 
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The derivatives dvi/dizj should be treated as 3 x 6 Jacobian matrices. The sys- 
tem ([9]) should be integrated simultaneously with the original equations ([1]) . The 
partial derivatives of v, are equal to the temporal derivatives of the corresponding 
partial derivatives of and can be obtained during the integration automatically. 
The initial values for dri/dirj and dvi/dTtj can be obtained by means of differ- 
entiation of the compound functions over the corresponding parameters. We 
do not give here the resulting expressions since they are rather complicated and 
numerous, and actually their main parts involving derivatives o f K r , K v over 7Tj 
can b e found in the classical literature on orbit refinement (e.g. iDuboshin et all 
Il976l sect. 3.3). The derivatives of ([7]) over m are also clearly calculatable, but 
they are rather complicated as well, and we have to omit them too. 



2.3 Numerical integrator 

We need some numerical integration method to solve the differential systems ([TJ 
and ([9|). We do not need it to be well-behaved (in any sense) on long time scales, 
since our integration timescale is rather short (~ 10 yrs or ~ 100 revolutions of 
the inner giant planet). But we do need it to be fast on this timescal e. After a 
few experim ents, we suggest that the integrators of the Everhart type ( EverharR 



Il973l 119741) would meet our requirements. We say "Everhart type" since we ac- 
tually used several integrators belonging to a general family of integrators similar 
to those constructed originally by Everhart. The original Everhart integrator was 
based on the Radau quadrature formula, b ased on certain as ymmetric splitting 
of the integration step. In accordance with lAvdvushevI |201oJj, it is possible to 
construct an integrator based on an arbitrary splitting of the integration step. 
These integrators represent, in fact, implicit Runge-Kutta integrators, equipped 
by an efficient method of evaluating the predictors on each step. Certain segment 
splittings (like the Radau one) allow to increase the order of the integrator sig- 
nificantly. Moreover, it is known that the use of the Legendre splitting makes the 
integrator symplectic. Symplectic integrators are very well suitable for long-term 
integrations, because they can preserve the Hamiltonian properties of the original 
iV-body problem over a much longer term (e.g., they show much slower energy 
error accumulation). It is important that in the case of the Everhart integrators 
the symplecticity can be obtained even with no significant trade-off or undesired 
side effects. 

The only significant disadvantage of such integrators is that they require con- 
stant time step for symplecticity. The algorithm itself does contain an optional 
step-size control mechanism, but variable step destroys the symplectic property. 
Therefore, we choose to use the same symplectic integrator both for the short and 
long integration terms. This is the 16-th order Everhart type integrator, based 
on the 8-node Legendre splitting. For short-term integrations (needed to perform 
orbital fits) we use the variable step. For long-term integrations (stability tests), 
we hold the step size fixed. 

The algorithm that we use is largely based on the Avdyushev's FORTRAN 
program cod(H with several modifications, aimed to increase its performance. We 



2 See preprint at http://www.scharmn.narod.ru/AVD/Gauss_15-2.pdf, in Russian 

3 See http://www.scharmn.narod.ru/AVD/Gauss-15.for and 



http://www.scharmn.narod.ru/AVD/GA USS- 32. for 
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omit the detailed description here, since we plan to provide it in a separate work 
in the future, along with the program source code. 



3 Radial velocity data for GJ876 and its preliminary planetary orbital 
solution 



The main datasets available are 162 Keck RV measurements from (Rivera et al 



l201Cft. and 52 HARPS RV measurements from (|Correia et alll2010|) . ICorreia et al.l 



(|2010T) also use RV datasets obtained at ELODIE and C ORALIE spectrographs, 
and there is also Lick RV dataset in ( Marcv et all l200lh . We include these data 



too, although it is necessary to note that they are considerably less accurate and in 
fact have insignificant effect on the system orbital solution — only a few per cent 
of the parameters uncertainties. The average internal RV precision of the HARPS 
and Keck datasets is about 1 and 2 m/s, their time span about 4.7 and 12.6 years, 
respectively. The average RV uncertainties of the Lick, ELODIE, and CORALIE 
datasets are about 10 m/s or more, and they span 6 and 8 years and 3 months, 
respectively. 

The internal errors stated in the published RV tables do not yet constitute the 
full RV uncertainty. It is well-known that high-precision RV exoplanet searches 
suffer from the phenomenon of the so-called RV "jitter" , which adds extra irregu- 
lar variations in their RV data. This excessive jitter was initially explained via star 
astrophysical a ctivity effects, but later it was shown that various systematic instru- 
mental effects ( Baluevl . 120091) may be significant as well (and may even dominate 



sometimes). Anyway, this jitter has to be taken into account during any further 
analysis. It is important here that in practice the effective value of the apparent 
RV jitter of the same star is usually different for different instruments, imply- 
ing different (and poorly known) weighting coefficients for different datasets. To 
take this jitte r into account properly, we utilize the maximum-likelihood approach 
suggested in (jBaluevl 120091) . This approach is based on maximizing the datasets' 
joint likelihood function, which depends on the RV curve (planetary) parameters 
as well as on the parameters determining the statistical structure of the errors 
in the RV data. This allows to estimate all necessary jitter values simultaneously 

with planetary pa rameter s. 

According to (jBaluevl [20083) . we will obtain the best fitting values of all in- 



volved parameters by means of maximizing the objective function 

ln£ = -Vvfln £ r full ,, i + ^-f-^- N ) ) - NbiVto, (10) 

where Nj stands for the number of the data points in the jth dataset, J de- 
notes the number of the datasets (J = 5 in our case), TV = Nj is the total 
number of observations, is the ith RV residual (O-C) in the jth dataset, and 
<7 fuii ji = a * j + °meas ji 1S the f un variance of the corresponding RV measurement, 
which incorporates the internal measurement variance c roeas as well as the jitter 
variance cr*j. The quantity 7 (not to be mixed with gammas in ([3}, which are 
indexed) is equal to 1 — d/N, where d is the number of free parameters of the RV 
curve (number of degrees of freedom in the RV model). This correcting divisor 
helps to remove the systematic bias in the RV jitter estimations, as discussed in 
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Fig. 1 The points mark the differences between the old llRivera et all [2005) and new 
( Rive ra et all [2010) Keck RV measurements for G J876, phased to the one year period (unit 
means end of a year). These differences show clear systematic variation of about 8 m/s in total, 
which is well above the typical measurement uncertainty of 3 — 4 m/s (old data) and 1 — 2 m/s 
(new data). This va riation was probably caused by errors in the data reduction pipeline used by 
iRivera et all j2005fl . Th e solid curve re presents the graph of the annual systematic variation, 
which was estimated in (Baluev, 2008b) on the basis of the old Keck data. 



This is necessary, because the residuals to the best fitting model 
are always systematically smaller than the original data errors, and without any 
correction they would yield underestimated value for the jitter. 

The function (|10p depends on the RV curve parameters (via rji, which involve 
the RV model) and on the RV jitter values (via o"f u njj). The values of these 
parameters, where the function (|10[) reaches its maximum, represent the necessary 
best fitting estimations. T o estimate a quality of a given fit, we use the statistic / 
defined in (|Baluevl,[2u09l) . This quantity represents a monotonous function of £, 
but is measured in m/s and thus is intuitively more clear and comparable to more 
traditional measures like r.m.s. (which we will use too, for some comparison). 

Other details of this method can be found in (|Baluevl . [2008cll2009t) . Here we 
use the same formalism, with the major difference that we now deal with the 
full gravitatio nal TV-body RV model, instead of the multi-Keplerian one. Also, on 
contrary with (Baluev, 2008c), we do not add in the RV model the terms describing 
possible annual systematic er rors in the da t a. Alt hough such annual errors existed 
in the old Keck data from ( Rivera et all 120051) , it seems the new data should 
have been fully corrected (see Fig. [l}. The HARPS observations for other stars 
seemingly were always free from such errors. 

For each planet, we will fit the following osculating orbital parameters: orbital 
period P (equiv. to the mean motion), mean argument of latitude tp, eccentric- 
ity e, pericenter argument u. We also fit common orbital inclination i (assuming 
coplanar system). In addition to the orbital parameters, we should also fit the 
planet masses. However, it is not convenient to adopt planet masses as primary fit 
parameters. In the traditional non-perturbed framework, the masses themselves 
are not determinable at all. In such a case, the RV oscillation semi-amplitude, 
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K, is fitted (as primary parameter), which allows to derive the minimum planet 
mass msini. It is well-known that without detectable interplanetary gravitational 
perturbations the inclination i remains unconstrained, and the true planet mass 
m remains unknown. In our case the system inclination is well-constrainable, but 
only for coplanar model. We will consider non-coplanar configurations below too, 
where individual inclinations may be poorly determined. In such a case it would be 
better not to mix the uncertainties of the RV amplitude and inclination, otherwise 
we risk to deal with significant troubles during numerical fitting, due to strong 
correlations between various fit parameters. Also, we usually will not include the 
innermost planet in the JV-body integration, assuming that on the observation 
timescale its motion is close to a Keplerian orbit. For this planet, we just have to 
leave with the RV semi-amplitude. 

Therefore, it would be better still to adopt the RV semi-amplitudes as primary 
fit parameters, instead of the planet masses. But then we must clarify what is 
"RV semi-amplitude" in the perturbed case. Actually, we cannot strictly define 
any "amplitude" for the perturbed motion, since this motion is not periodic. Nev- 
ertheless, in the case of GJ876 the deviations from the strict periodicity are small, 
and we still can determine the RV amplitude approximately, within the error of 
O(ni). Eventually, we need just to bind this RV amplitude parameter to the cor- 
responding planet mass. Therefore, we can just define it via mj using some simple 
formula close in shape to the formula from the Keplerian RV case. We adopt the 
following definition: 



K\f\ 



1/3 



^= sin* \W) ■ 

We must emphasize that this formula no longer expresses planet mass via the cor- 
responding RV semi-amplitude, as it would be in the case of unperturbed motion 
(/Ltj >C 1). Rather, it now defines the "RV semi-amplitude" via the planet mass. We 
may think of (|lip as of a definition of the osculating RV semi-amplitude, which 
completes the set of usual osculating orbital elements. So-defined "osculating RV 
semi-amplitude" is in fact just an intermediate fit parameter needed to separate the 
uncertainty of the orbital inclination from the uncertainty of the planet mass. The 
apparent semi-amplitude of the star's RV oscillation, caused by the corresponding 
planet, should be very close to the value of K defined in (|11[) . Furthermore, we de- 
cide to use not ev en the semi-am plitude K itself, but the quantity K = K\/l — e 2 , 



as it was done in (jBaluevl . [2008c) . Such choice allows to eliminate the eccentricity 
from the relation (|11[) and thus facilitates the conversions between fj, and K. 

To obtain some basic preliminary esti mations of th e syste m parameters, we 



take the Keck-only orbital solution from ( Rivera et all l2010h as a starting ap- 



proximation and perform the non-linear maximization of the likelihood function, 
as it was explained above. Table [T] contains the resulting best fitting estimations 
of all planetary parameters and RV jitter for different datasets. Note that here 
we exclude the innermost planet d from the integration, assuming that it moves 
along a Keplerian orbit in the common system orbital plane. Its orbital period is 
considerably smaller than periods of other planets, as well as its mass, and there- 
fore it does not show significant dynamical interaction with other planets on the 
observational timescale. Taking this planet into JV-body integration slows down 
the calculations dramatically. To ensure the gravitational influence of this planet 
is insignificant, we performed a similar (much longer) calculation, based on the full 
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Table 1 Best fitting coplanar orbital solution for GJ876 system (epoch JD2452000). 



parameter planet b (*) planet c (*) planet d planet e (*) 

fitted planetary parameters 



P [days] 


60.9904(68) 


30.1829(63) 


1.937886(18) 


124.51(52) 




K I'm /el 
l\ [ill / o\ 


213.21(34) 


84.65(36) 


6.18(29) 


3.41(33) 






341.13(20) 


71.09(46) 


357.6(3.2) 


299.3(7.3) 






0.0328(13) 


0.2498(28) 


0.178(44) 


0.008(27) 




I [°] 


248.7(2.9) 


252.08(51) 


224(16) 


181(77) 




* [°] 




56.1(1.5) 








derived planetary parameters 






m [M Jup ] 


2.377(42) 


0.747(13) 


0.0218(11) 


0.0482(47) 




a [AU] 


0.211018(16) 


0.131727(18) 


0.02110625(13) 


0.33961(94) 






RV data series and 


general fit parameters 






Keck 


Lick 


HARPS 


ELODIE 


CORALIE 


co [m/s] 


50.95(27) 


-31.4(4.9) 


-1337.87(42) 


-1864.1(3.7) 


-1904.0(4.6) 


o* [m/s] 


2.37(22) 


-11.3(6.0) 


1.63(23) 


21.3(3.2) 


19.7(4.3) 


r.m.s. [m/s] 


3.00 


27.6 


1.84 


33.5 


32.0 




1 = 


5.842/2.777 m/s 


, d = 26 







Each estimation is accompanied by its uncertainty in parenthesis (e.g., 0.30(10) means 0.30 ± 
0.10, and 30.0(1.0) means 30.0 ± 1.0). These u ncertainties were calculated from the Fisher 
matrix of the likelihood function {Baluev, 2009|). The uncertainties of the planet masses and 
semi-major axes do not incorporate stellar mass uncertainty. Planets included in the iV-body 
integration are marked with (*). The negative value for the Lick RV jitter means (symbolically) 
that the corresponding value of a\ is actually negative. The second value of I refers to the same 
fit based on only Keck and HARPS data (which offers practically the same estimations within 
a few per cent of the uncertainties). 



four-planet A?-body model. Almost all of the resulting best fitting parameters were 
practically identical, with negligible offsets of no more than 3% of the correspond- 
ing uncertainties. The only exception was the period of this same planet d, which 
decreased by 1.9 ■ 1CP 5 day, i.e. roughly by its uncertainty. From the statistical 
view point, such shift should not be neglected, since it would correspond to rather 
large one-sigma significance leveljj] This shift is nevertheless extremely small. It 
does not affect the motion of the other planets at a measurable level. Since the 
dynamical effects due to the innermost planet are so small, we neglect them below, 
and integrate only the system of three remaining planets. 

The apparent orbital period is also shifted due to the planetary aberration 
(also known as Roemer effect), which is caused by the fini t e li ght speed and is 
similar (in origin) to the Doppler effect ( Ferraz-M ello et all 12005)) . GJ876 radial 
velocity of -1.3 km/s (HARPS) implies that this shift should be -0.9 ■ 1CT 5 day 
(this is a "blue" shift, since the star is approaching). The cumulative correction to 
Prf, due to the perturbations and Roemer effect, is —1.0 ■ 10~ 5 day, which is about 
half of the corresponding statistical uncertainty. This correction should be added 



4 Throughout this paper, we will usually use the popular "n-sigma" style to specify various 
confidence probabilities. For example, we will say that a given statement is valid at an n-a 
significance level, when corresponding confidence probability is equal to the probability for a 
Gaussian random variable to deviate from its mean by no more than n times its standard de- 
viation ("sigma"). The significance levels of 1, 2, 3-<r correspond to the confidence probabilities 
of, respectively, 68.3%, 95.4%, and 99.73% (closer to 100% means more significant). 
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postfactum to all estimations of Pd in the fits in this paper, if such precision is 
necessary. 

So far, we limit ed ourselves to the coplanar four-planet system model used by 
iRivera et al.l (|20f0|) . When testing more complicated models, we found that the 
orbital fit from Table [1] shows statistically significant improvement if we add a 
free linear trend to the RV curve model. The estimated magnitude of this slope 
is about 0.18 ± 0.08 m/(s-yr), which results in a quite measurable RV offset of 
~ 2 m/s, accumulated during the observation time span. The formal sig nificance 
of th is trend, as derived from the corresponding likelihood ratio statistic ( Baluevl . 
l2009h . is about 2.3a. Even though we have not yet took into account several im- 
portant effects that we will discuss in subsequent sections, such significance is too 
large to be neglected without investigation. Initially, we interpreted this long-term 
slope as the geometrical secular RV acceleration effect, mentioned bv lCorreia et al.l 
(2010), which should be equal to 0.15 m/(s-yr) for GJ876. In fact, subtracting this 
predicted slope from the RV data allows to get rid of any significant trend in the 
fitted model. Actually, we started our research based on such corrected data, no 
longer bothering about any RV trends. However, A. Correia and E. Rivera later 
confirmed (in private communication) that both published RV time series already 
have this trend subtracted off. Therefore, we need to find another interpretation. 
We can see three equally plausible sources: another long-period unseen compan- 
ion of the star, a tiny long-term instrumental drift, or some extra errors in the 
RV reduction pipeline. Regardless of the actual nature of this probable slope, we 
should try to take it into account, since it could significantly affect our results. 
Since its source and magnitude remains a priori unclear, we add a free linear term 
to the RV model. Most of our results below will refer to this model with trend, 
although sometimes we will make a comparison with the trend-free model. We 
will also return to a more rigorous estimation of the actual significance of this RV 
trend in further sections. 



4 Correlated radial velocity jitter 

Let us first investigate whether the currently available RV data show some residual 
periodicities, in addition to the current four-body model of the planetary system 
and possible linear trend. To do t his, we utilize the common periodogram-based 



approach with modifications from (|Baluevl . 120091) . The main modification is neces- 



sary to take into account the likelihood function (|10[) , whereas the traditi onal pe- 



riodograms (e.g. the Lomb-Scargle one) implicitly utilize the \ 2 function ( Baluevl 
2008a). Each value of the periodogram that we are about to use, is associated with 
the likelihood ratio statistic, measuring how much our RV model improves, when 
we add a probe sinusoidal signal to it. Another im portant modifi cation, which we 
must highlight here, is that the periodogram from ( Baluevl l2009h is not the tradi- 



tionally used periodogram of the fixed RV residuals to the base best fitting model. 
Rather, it requires a full re-fitting of the whole set of the parameters and full 
re-evaluation of the RV residuals for each periodogram value. This increases the 
sensitivity of the periodogram to faint periodicities, since we can better model the 
cumulative RV variation, when the probe signal indeed exists. We will refer to such 
periodogram as "residual periodogram" , on contrary to the "periodogram of the 
residuals" . The residual periodograms are obviously much more computationally- 
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Fig. 2 Raw and smoothed residual periodograms for the four-planet model of GJ876 (with 
eccentricity e e always fixed at zero) plus free linear RV trend. Top and bottom plots show 
periodograms for the HARPS and Keck datasets. The smoothing was performed using moving 
average over the frequency segments of 0.09 day -1 . The periodograms (in the Keck case 
especially) show excessive power at the frequencies / < 0.2 day _1 and / > 0.8 day -1 and a 
relative depression in the middle of the segment. 



demanding than the traditional ones, especially when dealing with Newtonian 
iV-body fits. 

The joint residual periodogram of all available RV data, corresponding to the 
base model from Table [T] shows no isolated peaks above the apparent noise level. 
However, it does not look like an usual white noise periodogram as well. On con- 
trary, the noise level itself demonstrates clearly varying structure. This variability 
becomes especially clear, when we plot the periodogram in the linear frequency 
scale, instead of the logarithmic one. Fig. [2] shows such periodograms, plotted 
separately for the Keck and HARPS datasets (i.e., assuming the probe signal is 
present in only Keck or only HARPS data). The joint periodogram represents some 
mixture of these. We can see that in case of the Keck data there is an excessive 
power at low frequencies (long periods) and near the unit frequency (period close 
to one day), with a depression in the middle of the segment. A similar frequency 
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distribution, although somewhat obscured by irregular variations, can be seen in 
the HARPS data too. These frequency spectra are very resistant with respect to 
various modifications in the RV curve model. The plots from Fig. [2] correspond 
to the circular orbit of the fourth planet, but assuming other reasonable values of 
e e and uj e (see the next section) does not significantly affect the smoothed peri- 
odograms (although individual periodogram peaks can be affected). The influence 
of the long-term RV trend on these spectra looks also insignificant, as well as the 
influence of possible system non-coplanarity. We could not find a way to explain the 
non-uniform shape of the smoothed periodograms via any possible shortcomings 
in the RV curve model. Therefore, we need to consider another explanations. 

The errors in astronomical time series are usually assumed mutually uncorre- 
cted. Such uncorrelated sequence of errors is also known as white noise, which is 
called so because of its uniform frequency spectrum. A non-uniform spectrum in- 
dicates autocorrelated residuals, according to the Wiener-Khinchin theorem. This 
means, in particular, that the fitting methods that we used above, actually are 
not applicable here, since they all are based on the assumption of uncorrelated 
RV errors. In such a case, the correlated RV noise could be misinterpreted as 
some deterministic variation, and could cause unqualified systematic errors of the 
estimations in Table [T] 

The variations of the noise level in Fig.[2]are rather large. The max/min ratios 
for the smoothed periodograms are 3.1 and 7.4 for the Keck and HARPS peri- 
odograms, respectively. Nevertheless, the apparent variations of the periodogram 
noise could also emerge purely by chance, because of the finite number of observa- 
tions (even when the original noise is white). To demonstrate that the variations 
of the average noise level in Fig. [2] are statistically significant indeed, and to check 
in which fraction they are significant, we carry out some Monte Carlo simulations. 
We adopt the fit from Table [1] as the basic one, and carry out a series of boot- 
strap simulations of the residual periodograms plotted in Fig. [2j The bootstrap ran- 
dom shuffling destroys any possible correlations between different residuals, so the 
simulated RV noise is practically white. Therefore, the simulated periodograms 
should contain only natural random variations, which we should expect for the 
uncorrelated error noise of the same variance and distribution. Each simulated pe- 
riodogram was evaluated using literally the same algorithm as real periodograms 
in Fig. [2l and was further smoothed to access its average level. The smoothing was 
done using moving average over frequency segments of w 0.09 day" 1 . The result 
of this procedure is a bunch of simulated smoothed white noise periodograms. 
Then we calculate confidence ranges for the smoothed periodograms, based on the 
simulated periodogram set (separate confidence range for each frequency). Also, 
we calculate for each smoothed simulated periodogram the ratio of its maximum 
value over the whole available frequency range to the minimum one, and count how 
frequently it exceeds the same ratio of the original periodogram of the real data. 
It will help us to check, whether the observed variations in the basic smoothed 
periodograms are statistically significant or not. 

The results of these simulations are plotted in Fig. [3j We can see that the 
simulated range of the smoothed periodograms remains almost constant in fre- 
quency. This simulated range is quite narrow for the Keck case, and more wide 
for the HARPS case, obviously because the number of the Keck observations con- 
siderably exceeds the number of the HARPS ones. Both periodograms of the real 
data do not stay in the simulated confidence ranges. In the Keck case, none of 100 
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Fig. 3 Smoothed residual periodograms for the four-planet model of GJ876 in comparison 
with their simulated levels expected for uncorrclatcd noise. Top and bottom plots correspond 
to the HARPS and Keck datasets. The smoothing was performed using moving average over 
the frequency segments of 0.09 day -1 . The graphs show the actual smoothed periodograms, 
their simulated mean levels and levels corresponding to the 16% and 84% percentiles (two-sided 
one-sigma limits). The max/min ratios of the original smoothed periodograms are also printed 
in each panel, along with their simulated one- and two-sigma upper limits (in parenthesis). 
We can see that both Keck and HARPS data show statistically significant deviations from the 
white noise in terms of their power spectra. 



simulated smoothed periodograms could demonstrate the same or larger max/min 
ratio as we see in the corresponding real periodogram. This means that the noise 
in the Keck data is not white, with high confidence probability well above 99%. In 
the HARPS case, the confidence probability is smaller - approximately 95% - but 
nevertheless is high enough to say that a similar non-whiteness probably exists 
in the HARPS data too, although significantly obscured by the normal random 
variations. 

We will not try to determine here possible physical sources for such behavior 
of the RV noise, since this is an topic for another research. Regardless the actual 
sources, two main practical questions arise now: how much the noise non-whiteness 
affects the best fitting orbital configuration from Table [TJ and how to correct this 
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effect? These problems must be solved before we can go any further. We find 
that the correlated noise i s already a routine issue in the exoplanetary transit 
searches ( Pont et al.l 12006'). The photometric observations in these surveys often 
contain the so-called "red" noise component. Such name is due its power spectrum, 
which monotonously decreases with the frequency, thus making long-term ( "red" ) 
variations to prevail over the short-term ("blue") onet[f|. Actually, the red noise 
is rather common phenomenon: it is not limited to only astronomy, and it is 
frequently faced in very different branches of the science. 

The red noise indeed can easily explain the shape of the periodograms in Fig. [2) 
the low-frequency "hill" is just directly seen in both graphs, and another (smaller) 
hill near the unit frequency is its alias (caused by severe diu r nal g aps in the 
time series). The transit observations considered in ( Pont et al.l . I2006T) allowed a 
simplified red noise reduction, due to their very specific distribution in tim e and 
a specific of their photometric transit models. In fact, IPont et"al ] (|2006f) could 
avoid, for instance, any assumptions about th e shape of the nois e autocorrelation 
function. Unfortunately, the approach used bv lPont et al.l ( 20061) is not applicable 
to our case. We have to model the correlated noise in full. 

The details of the algorithm, which we use to eliminate the impact of the red 
RV jitter, are given in Appendix [A"l Here we only note that this algorithm is based 
on a certain modification of the function (|10[) . to take the correlated noise into 
account. After the maximization of this new objective function, we get (at once) 
the estimations of the RV curve parameters (therefore, of the planetary orbital 
elements and masses), and several RV noise parameters: different "white" jitters 
for separate time series, er w hitej, the common "red" jitter shared among different 
time series, a le d, and the red noise correlation timescale r. Such noise separation 
could be realistic if the red noise was actually caused by the star (e.g. by some 
long-living photospheric phenomena) rather than by the instruments. 



5 Determinability of the fourth planet orbit 



Further investigation of the likelihood function (|10|) showed that it possesses mul- 
tiple maxima, which do not differ much from each other, in terms of the goodness- 
of-fit measures like r.m.s. or the statistic I. Such phenomenon could indicate that 
the whole orbital model is actually ill-determined due to the lack of the RV data 
and/or its in sufficient precis ion, like that was for the extrasoalar planetary system 
discussed in (|Baluevl . l2008dY However, it is not the case here. For the GJ876 plan- 



etary system, the irregularities of the likelihood function are all related to only two 
of the free parameters from Table [1] namely the eccentricity and the pericenter 
argument of the fourth planet. Fixing them at some reasonable a priori defined 
values makes the shape of the likelihood function very regular and practically 
parabolic, as it should be when the RV model is well-linearizable with respect to 
the remaining free parameters. 

We check this for the following pairs of parameters: (e c ,i), (P;,,Pc), (i/^, V'e), 
and (lo{,,ujc)- The estimations of these variables possess the largest mutual cor- 
relations in the pairs, respectively —0.98, —0.82, —0.81, and —0.72 (for e e fixed 



5 There is also a narrow notion of "red noise" as a synonym of the Brownian noise, that has 
a specific frequency spectrum <x l// 2 . We understand the term "red noise" in a general sense, 
however. 
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at zero). Correlated pa rameters shoul d be usually more affected by non-linearity 
effects, as discussed in (lBaluevll2008cV For each above pair of parameters, we plot 
two-dimensional contours of the function (| lUj) . We choose three contours which 
outline confidence regions for the selected parameter pairs, corresponding to the 
asymptotic (N — > oo) confidence probabilities of l,2,3cr. The necessary confidence 
probabilities were calculated from the modified likelihood ratio statistic (logarithm 
of the differenc e between the maximum and a given value of (|T0j| ) , as discussed in 
(|Bamevl l2009). After that, we compare these confidence regions with the results 
of numerical Mo nte Carlo simulation s. For this goal, we use the bootstrap method 
as described by iMarcv et all (|2005l) . In this method the simulated "errors" are 
generated by means of the random shuffling of the RV residuals to the best fitting 
RV curve model. This method of Monte Carlo simulation is rather widely used, 
because it does not require any assumptions about the shape of the error distri- 
bution, and is expected to work even for non-Gaussian errors. We introduce only 
two relatively cosmetic modifications to this method. First, during the shuffling 
we do not mess the residuals belonging to different datasets (we shuffle different 
datasets separately from each other). Second, after each shuffling trial, we rescale 
the residuals according to the total variances for the corresponding observations 
(which include RV jitter). We consider only the uncorrelated noise model, since 
random shuffling destroys any correlations anyway. 

During all of the activities described in the previous paragraph, we always 
held e e and u e fixed at zero, to eliminate all non-linearity effects associated with 
these parameters. The results of these calculations are shown in Fig. [3] We can 
see that all of the resulting confidence domains have almost elliptic shape, and the 
simulated sets of points are always in very good agreement with these confidence 
regions. This means that when the parameters e e and ui e are fixed, the remaining 
parameters behave almost as in the linear least-squares problem: there is a single 
maximum of the likelihood function (|10[) , and this function has almost parabolic 
shape in the vicinity around the maximum. The estimations of the parameters 
possess almost Gaussian distribution, and their uncertainty estimations are pretty 
reliable. Notably, even the pair involving ip e shows no significant non-linearity 
effects, despite rather large uncertainty in this parameter. 

The behavior of the pair (e e ,oJe) is severely different. The similar confidence 
contours for these variables look very irregular (Fig. [5j) and outline two main local 
minima of I. The first local solution is the one listed in Table [T] and the second 
one has larger eccentricity e e ~ 0.12. This second solution appears actually a bit 
more likely, when we use the white noise model. But for the model with correlated 
noise, the first solution becomes the leading one. The picture is also sensitive to 
the trend term in the RV curve model. Therefore, this multi-extrema structure is 
model-dependent. Can we trust to this fine structure of the likelihood function at 
all? Most probably, we cannot: it may easily change even further, as some other 
non-traditional RV noise effects are discovered. The irregular details in Fig. [5] are 
driven by RV variations at the level of ~ 10 — 30 cm/s, which is well below the 
internal measurement errors (~ 1 m/s). Such small variations are hardly related to 
the actual radial velocity of the star. They are more likely caused by unqualified 
systematic errors or systematic astrophysical noise in the data, which could easily 
have amplitude about 1/3 of the random measurement noise. Therefore, the only 
reliable information that we can obtain here is that the eccentricity e e probably 
does not exceed ~ 0.15 — 0.20, and uj e values near ~ 45° are more favored than 
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Fig. 4 Predicted confidence regions for the several selected parameter pairs of the GJ876 
planets, in comparison with numerical simulations. Each panel shows three contours of the 
function HlOt , which were calculated to outline the confidence regions corresponding to 1, 2, 3- 
ct significance levels, in the framework of the linear least-square problem (i.e., asymptotically 
for N — > oo). These regions are all almost elliptic and are in good agreement with results of 
the bootstrap simulations, shown as points. This indicates that the problem is well lincarizablc 
with respect to the selected parameters. The eccentricity e e was always fixed at zero here, to 
eliminate the non-linearity effects coming from the (e e ,u} e ) pair. Note that these plots only 
demonstrate the linearity of the selected parameters, but they should not be considered as a 
source of information about the corresponding uncertainties, since many effects are not taken 
into account here yet. 



the opposite ones. All other fancy details in Fig. \5\ represent just some misleading 
noise component, which we must try to eliminate. 

This pseudo-detailed structure is not caused by some specific property of the 
fitting method applied here. Some points in Fig. [5] indeed yield smaller scatter 
of the RV residuals than the others, and this pattern is inevitably irregular. Any 
other classic fitting method that use the function (|10[) or any similar function as 
a single goodness-of-fit measure, would yield similar results. To understand the 
source of such irregular behavior of only two of the parameters, we need to trace 
where the information about these bad-behaving parameters comes from. To do 
this, we first replotted Fig. [5] for the case when the planet e is excluded from the 
iV-body RV model, and its contribution was modeled by a Keplerian function. Such 
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Fig. 5 Likelihood function contours for the parameters (ecosaj, esincj) of the planet GJ876 e, 
assuming the RV noise is non-correlated (white) and correlated (white+red). The gray points 
along radial lines mark the orbital configurations disintegrating in less than 10000 yrs. Each 
panel shows four contours of the likelihood function, which were calculated to outline the 
regions of the asymptotic confidence probability matching the 1, 2, 3, 4-tr significance levels. 
However, the complicated irregular and model-dependent structure of these regions indicates 
that they are unreliable and may be driven by extra unqualified systematic errors in the RV 
time series (see text for discussion). 



plot was already pretty regular and elliptic, just as those in Fig. [3] However, the 
corresponding confidence domains appeared much wider than in Fig. [3] allowing e e 
values of up to 0.3 — 0.4. This means that the contours in Fig.[5]could not be settled 
by the non-sinusoidal shape of the RV variation contributed by the fourth planet, 
as it would be in the non-perturbed Keplerian case. Instead, the information about 
the likely values of e e and cu e is extracted mainly from the perturbational effects, 
which the fourth planet imposes to the motion of other planets. Therefore, the 
main source of the irregularities in Fig. [5] is the dynamical interaction of the 
planet e with other planets in the system. This dynamical interaction allows to 
considerably shrink the region of admissible values of e e and uj e , but by the cost 
of considerably irregular dependence of the RV curve on these parameters. As a 
result, we see some unreliable structures inside this region. 

In view of the unreliable behavior of the parameters e e and uj e , the best course 
of action for us will be to recognize that all points within some wide enough contour 
in Fig. [5] are equally likely. To determine, which contour should serve as a bound- 
ary, we apply a dynamical stability test to all orbital solutions spanning Fig. [5] 
Each point in this graph corresponds to some best fitting orbital configuration of 
the system. For each of these configurations, we perform the numerical integration 
over 10 4 yrs and check, whether a given configuration disintegrates during this 
test term. The integration term of 10 4 yrs is actually pretty short in comparison 
with cosmological timescales, but nevertheless it contains about thousand of sec- 
ular periods of the system. Therefore, it roughly corresponds to ~ 10 7 yrs, when 
rescaled to the subsystem of the giant planets in our Solar System, for instance. 
Such time segment is long enough to determine approximate stability boundaries. 
In Fig. [5l this stability boundary passes close to the formal Aa confidence contour 
(in case of the white+red noise model) . Since the stability region should somewhat 
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Fig. 6 Residuals of the best fitting three-planet model of GJ876, phased to the orbital period 
of the fourth planet. Small crosses mark the Keck data, fat points stand for the HARPS data. 
Solid curve is the unperturbed sinusoidal model of the RV oscillation due to the planet e. 



shrink, when the integration time increases, we believe that it is safe to adopt the 
formal 3a contour as a boundary outlining the admissible values of e e and w e - 

In view of such rather poor determinability of the fourth planet eccent r icity, one 
can ask: whether this planet exists at all? It was detected b y lRivera et al.l (|201fj|) on 
the basis of the Keck data only, and ICorreia et all ( 2010^) did not note it in their 
HARPS and old Keck data. We also checked the residual periodograms for the 
three-planet model, and find that the fourth planet is indeed strongly supported 
by the Keck data, but the HARPS periodogram does not show a significant peak 
above the apparent noise level. One can suspect that the RV signal from this fourth 
planet could be actually caused by some spurious drifts in the Keck data. 

However, these doubts dissolve when we look at the RV residuals to the three- 
planet model directly. We can see (Fig. [6]), that the available HARPS data actually 
confirm the RV signal from the fourth planet and are in good agreement with the 
Keck data. This agreement just is not statistically significant yet, however it may 
become more significant, when more HARPS data are accumulated. At present, 
we have no observational basis for doubts in the existence of the planet GJ876 e. 



6 The reference orbital fit 

Using the algorithm from Appendix [A] we can now obtain the full set of the 
parameters, taking into account the effect of the correlated noise. However, we 
must also address the influence of the bad-behaving parameters (e e ,uj e ). Without 
any extra care, we will get just a few similar orbital solutions having unrealistic 
uncertainties. We adopt the following approach. Since other parameters behave 
almost linearly (see Fig.U]), we first perform a basic orbital fit, fixing the value of e e 
and u! e at some reference realistic value (say, e e = 0) . The resulting error estimation 
for each of the fitted parameters reflects only a fraction of the full uncertainty. We 
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Table 2 Reference orbital solution for the GJ876 system, assuming only white RV noise 
(epoch JD2452000). 



parameter planet b (*) 



P [days] 


60.990 | 


'e, ; 


H20A 
-30 J 


K [m/s] 
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33, 


+21 
-29 


V> [°] 


341.11 ( 


19, 


+5l' 
-62 
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0.0332 | 
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+ 17' 
-38 


w [°] 
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+9.3' 
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fitted planetary parameters 
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84.48 (35, ',' ! (>.20 j 2> 
72.0 (o.4, T i i ) 357.5(3.1 



+9 

-16 
+ 76 
-61 
f 1.1 
-0.7 



1.937888 1 18, +°) 



+14 
21. 



HO. 5 
-0.5 , 



0.2511 (28, +42) 

+o.e' 
-1.6 , 



0.148 (45, 
217 f 19,^1) 



+1 
-3.8 



252.2 (0.5, 

56.1 (l.5, 
derived planetary parameters 
0.750 (13,+^) 0.0221 (ll, 
0.211021 (l5, +11) 0.131726 fl7, +H) 0.02110627 fl3, +? 



2.39 4. 



+12 



planet e (*) 



5.4 

6, 



124.5 (0.4, 
3.6 (0.3, 
300 (6, +£) 
0(< 0.19) 
0(unconstr.) 



0.051 5 



+4 1 
-14, 



0.3397 ( 7, 



Keck HARPS 
RV data series and general fit parameters 



c [m/s] 

ci [m/(s-yr)] 

o-whito [m/s] 
r.m.s. [m/s] 



50.66 28, 



-1338.71 53 



0.173 I 74, +£) 



1.61 22, 



+60 
-80 



—19 



2.31 I 22, 
2.97 1.81 
= 5.780/2.764 m/s, d = 27 



— 16 
-18 , 



The same notes as in Table [T] apply here, except for a more complicated treatment of the 
parameter uncertainties. Each estimation is now accompanied by three uncertainty values in the 
parenthesis: the first value denotes the la uncertainty due to the usual linear statistical effects, 
assuming e e fixed at zero. These values were calculated in the traditional way. The remaining 
pair of values in the parenthesis (one above another) reflect the asymmetric uncertainty inferred 
by the non-linear parameters e e and ui e . The last figure given in each uncertainty value always 
maps to the last figure in the corresponding estimation. The estimations themselves, as well 
as the values of the r.m.s. and I are given for the basic solution e e = 0. We omit the values 
related to the Lick, ELODIE, and CORALIE time series. See text for further details. 



also need to estimate the uncertainty caused by the non-determinability of the 
parameters e e and ui e . To do this, we vary e e and ui e inside their admissible region 
(which was defined in Sect. [5]) and see, how much this affects the best fitting values 
of other (well-behaving) parameters. This gives us the remaining part of the total 
parameter uncertainty (generally asymmetric). 

The results of these calculations are given in Table [2] (white noise model) and 
Table[3](white+red noise model). Note that when we varied e e and uj e , this mainly 
affected the estimations of other parameters themselves, and the estimations of 
the statistical uncertainties remained almost constant. The first (probabilistic) 
uncertainties for the most of the parameters in Tables [2l3l should be quite reliable, 
with a very few exceptions though. These exceptions are caused, however, by the 
statistically unsuitable non-linear parametrization of the planetary model, rather 
than by the internal non-linearity of the problem itself. The estimations of bounded 
parameters, like the eccentricity or the RV jitter (both > 0) are considerably non- 
Gaussian and asymmetric, if the formal uncertainty range can cover the forbidden 
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Table 3 Reference orbital solution for the GJ876 system, with red RV noise taken into account 
(epoch JD2452000). 



parameter planet b (*) 



P [days] 
K [m/s] 
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m [M Jup ] 
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Keck HARPS 
RV data series and general fit parameters 

c [m/s] 50.79 (33,+^) -1338.61 f 63, 

d [m/(s-yr)] 

o-whito [m/s] 
cr rcd [m/s] 

T r cd [day] 
r.m.s. [m/s] 



+45 
-83 



0.49 (54, +_ 3 9 4 ; 

+37\ 
16 J 



The same as in Table [2] but assuming that the RV noise contains a common red component 
shared between different time scries. 



values. It is nonetheless very easy to find a better (more linear) parametrization 
in such cases. In particular, it is better to consider the pair (e^cosct^e^sinc^) 
instead of (e^a^), as it will be done below. When the value of some RV jitter is 
close to zero (in comparison with its uncertainty) then it is better to deal with 
the squared ji tter instead, with the uncertainty properly rescaled (see Table 1 by 
iBaluevI (|2009|) ). 

Looking at these orbital fits, we note, at first, that the linear RV trend appears 
more disputable than it seemed before: some values of e e and uj e infer too low 
significance for ci (about 1 — 1.5a). Basically, the need for this trend can be elimi- 
nated by means of choosing appropriate values for e e and uj e (namely, in the right 
and top regions in Fig. [5]). On contrary, it is too early to claim that this trend does 
not exist at all. Other possible values of e e still require this RV trend. Allowing 
non-coplanar configurations (see Sect. [7]) keeps this RV trend almost intact, so it is 
not easy to explain this RV trend via orbital non-coplanarity as well. We conclude 
that this issue can be resolved by future observations. 

Second, the red noise correlation timescale r rec j is in fact poorly constrained. 
We can only say that r re( j has an order of days. On contrary, the magnitude of the 
red jitter itself, cr Tedl is rather well constrained and is well separated from zero. 
This suggests that although the very existence of the red noise in the data looks 
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Fig. 7 The confidence regions for the eccentric parameters of the planet d, plotted in the 
same way as in Fig. [5] The stability test integrations were not carried out here, because 
these parameters are not crucial for the system stability. The confidence contours (thick lines) 
were plotted for the 1,2,3-cr probabilities and assuming e e is fixed at zero. The remaining 
uncertainties, inferred by the bad-behaving parameters e e and ui e , are rendered as gray domains 
around the best fitting points for e e = (marked as crosses). See text for discussion. 



supported, the characteristics of this jitter are still difficult to assess. Neverthe- 
less, at present even a rough estimation of r re( j can provide important physical 
information. 

Third, we can see some drop in the eccentricity e d , after we take the red noise 
into account. Since the apparently eccentric orbit of the planet d has an important 
value for, e.g., the star-planet tidal interaction theory (e.g. iFerraz-Mello et all 
120081) . we should investigate the parameters e d and u> d more closely. Fig.[7Jcontains 
the confidence contours for the parameters (e^cosai^e^sinu;,;;), constructed for e e 
fixed at zero, like in Fig. 2] We can see that they are perfectly elliptic and, in the 
white noise case, also agree with the bootstrap simulations (not shown). 

Before drawing any conclusions, we need to characterize the uncertainty coming 
from the parameters e e and u e . Again, we adopt the 3a contour in Fig. [5] as a safe 
region of admissible values for e e and u> e , not paying attention to the apparent 
concentrations inside this region. Each point in this region refers to some orbital 
fit with fixed e e and ui e - The values of e d and uj d from these fits would mark the 
centers (best fit points) of the error ellipses in the plane (e^ cosw^, sinw^), if 
we actually constructed such confidence ellipses for each admissible (e e ,oj e ) value. 
For instance, two best fit points, shown as a crosses in Fig. correspond to 
the central points e e = in Fig. [5] Furthermore, instead of this single point in 
each panel, we can construct the full set of points (e^ cosu;^, sincj^), which are 
mapped from all admissible values of e e ,oj e - These sets are rendered in Fig. [7] as 
gray domains around the basic best fit points. These domains reflect how much 
the centers of the error ellipses in Fig. [7J may shift, while the values of e e and 
Lj e are varied inside the admissible region. Although these centers may shift, it 
turns out that the shape and size of the elliptic contours in Fig. remain fairly 
constant for different e e and to e . The picture only shifts in the fairly solid way. 
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This means that the uncertainties inferred by the bad-behaving (i.e., non- linear) 
parameters (e e ,to e ) are well-separable from the uncertainties inferred by other 
(well-linearizable) parameters. 

Therefore, we may treat now that the resulting uncertainty region in the plane 
of {e d cosoj^, e d sino;^) is constituted in a cumulative manner from the usual prob- 
abilistic confidence domains, inferred by the linear estimation theory with e e fixed 
at zero (shown as regular elliptic contours), and from the uncertainty domain in- 
ferred by the non-linear parameters e e and ui e (shown as less regular small gray 
regions). We consider that all admissible values of e e and cj e are equally possible, 
and therefore the gray uncertainty regions in Fig. [7] should be understood as struc- 
tureless solid entities, which just indicate how the original probabilistic confidence 
regions should be bloated to take into account the uncertainty coming from e e , w e . 

We can see from Fig. [7] that e d is indeed inconsistent with zero, when we 
analyze RV data assuming traditionally that the RV noise is white. The significance 
of this non-zero e d is well above the 2a level, even when we take into account 
the uncertainty inferred by the non-linear parameters e e and cj e . However, it is 
now obvious that this apparently significant value of e d is likely a result of the 
misinterpreted red RV noise. When the red noise is taken into account, the best 
fitting value of e d moves closer to zero, and simultaneously the uncertainty regions 
expand. Taking into account the uncertainty coming from e e and u e , we realize 
that e d is in fact consistent with zero at the significance level of hardly above 
la. The things remain similar when our RV curve model does not contain the 
linear RV trend. In that case, e d is non-zero at > 3a level for the white noise 
model, but for the correlated model this significance drops to the same ~ la level. 
Although we still cannot retract the values of e d as large as ~ 0.15, we nevertheless 
have no observational evidences that e d is actually non-zero. Further observations 
can eventually solve this question for sure, but at present it is too early to claim 
that the non-zero eccentricity e d was confirmed. The apparent non-zero value of 
e d reported in previous works represents basically an effect of misinterpreted red 
noise in the RV data. 



7 System non-coplanarity 

One of the primary goals of this paper was to characterize the mutual non- 
coplanarity between planets b and c or at least to put some limit on it. In the 
non-coplanar case, we have two separate variables for the osculating inclinations 
if, and i c , and also a pair of extra parameters determining the orientation of their 
ascending nodes, Qf, and Q c - Two latter parameters, however, cannot be esti- 
mated independently, because the problem is invariable with respect to arbitrary 
rotation around the line of sight. In fact, we can determine only the difference 
AQjj C = fi c — flfr. The main quantity that we are interested in, in view of the orbit 
non-coplanarity, is the mutual orbital inclination I for the planets b and c, which 
is a function of ib, i c , and Ail^c- 

We must note that during this work we will deal, most probably, with small val- 
ues of I, comparable to its statistical uncertainty. Such situation reveals many prac- 
tical difficulties, because our estimations of I will have significantly non-Gaussian 
distribution. Similar troubles arise for planetary eccentricities, when they are small 
enough (comparable to their uncertainties). In the both situations, the issue is 
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to observer 



Fig. 8 Illustration of the angle <P definition. For simplicity, both orbits are assumed circular 
and of the same radius. The direction Q is determined by the orbits intersection point where 
the second planet ascends over the first orbital plane (and not the opposite intersection point). 
The directions i?i and Q2 arc defined in the similar way: each planet should ascend over the 
sky plane in the point where the inclination angle is determined. Such definitions allow for all 
inclinations ii, 12, 1 to be always non- negative (keeping implicitly their signs in the orientation 
angles i?2,^)- See text for the detailed discussion. 



rather formal, however. The non-Gaussian behavior of such estimations is caused 
by the bad choice of the parametrization, rather than is inferred by the RV model 
itself. In case of the eccentricity, the problem is usually eliminated if we consider, 
instead of the eccentricity, the pair of parameters (e coscj, e sinw). Estimations of 
these parameters usually have almost Gaussian bivariate distribution, even when e 
is small. All apparent problems in this case are caused by the trivial singularity of 
the polar coordinate system (e,w) in the point e = 0. It is very likely that the mu- 
tual inclination 7 should be affected by a similar singularity at 7 = 0, which could 
be eliminated by means of transition to the variables like (sin 7 cos <P , sin 7 sin $) , 
where $ is an extra auxiliary variable, determining some extra orientation angle 
related to the mutual orbital inclination. Note that the variables 7 and <P should be 
independent in the sense that for any pair of 7 e [0, it] and <P € [0, 2tt] there should 
exist a valid orbital configuration. Due to that requirement, we cannot choose <P 
to be equal to, e.g., the angle between the planetary ascending nodes AO. For 
any I < n/2, the value of AS7 cannot exceed 7, and such non-constant constraint 
will imply nothing except extra difficulties. We could use some orientation an- 
gle in the Laplace plane of the system. Such choice would be physically justified 
and also symmetric with respect to the two planets involved, but it would mess 
the geometric parameters (like inclinations) with planetary masses, which is not 
desirable. 

We adopt the following purely geometric definition of the auxiliary angle <P. 
Given two abstract "first" and "second" planets, we can determine two mutual 
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orbital nodes of these planets. We choose the node in which the second planet 
ascends over the orbital plane of the first planet, and define $ as the orientation 
angle of this node in the plane of the first orbit, counting it from the usual ascend- 
ing node of the first orbit. This definition is schematically illustrated in Fig. [8j 
With a help of the classical spherical trigonometry, it is not hard to derive the 
formulae expressing the new osculating angles I,<P via the original parameters 
ii,i2, AO = O2 — 0\. They look like: 

sin I cos = — sinii cos 12 + cosii sin 12 cos AO, 
sin/ sin <P = sin 12 svnAO, 

cos / = cosii cos 12 + sinii sin 12 cos AO. (12) 

The inverse transition is given by the equalities: 

sin %2 cos Afl = sin i\ cos / + cos i\ sin I cos <£, 
sini2 sinZlJ? = sin I sin <P, 

cos %2 = cosii cos/ — sinii sin/cos#. (13) 

Obviously, it is impossible to express all three original orientation parameters 
via only / and $. We must also know i\ to obtain 12 and AO. Therefore, such 
definition of <P is not symmetric with respect to the two planets: the first orbit 
serves as a reference plane. In case of GJ876, we choose the planet b to be this 
"first" planet, since its RV amplitude is considerably larger, and thus its orbital 
plane orientation is determined with better precision. The planet c will be the 
"second" planet, and the orbit of the planet e is assumed to lie in the common 
Laplace plane of the system (at the epoch of osculation) . The planet d was also 
assumed to move in the Laplace plane, though this assumption only affects the 
mass estimation of this planet Q 

Given the formulae (|12I13[) . we can easily carry out constrained RV fits with 
/ and $ fixed at any desired values. Technically, such constrained fitting may 
be done, for instance, by means of expressing 12 and AO via /, <E>, and i\ (thus 
eliminating 12 and AO from the set of free parameters). Therefore, we can perform 
a series of such constrained fits on some regular grid of / and $ and then to plot the 
resulting likelihood contours in the plane, e.g. , (/ sin , / cos $) . Such plot basically 
visualizes the confidence regions for these variables, similar to those regions that 
we have already constructed in Fig. 2] 

Results of these calculations are shown in Fig. [9] for both noise models (white 
and white+red). It is notable that in the white noise case the inclination / shows 
some non-zero value at the significance level of 1.9a, although this significance 
becomes somewhat smaller when the uncertainty in (e e ,oj e ) is included. These non- 
coplanarity signs are further softened, when the red noise is taken into account. 
In this case, / is consistent with zero at 1.2cr level, even when the e e and co e 
uncertainties are neglected. Eventually, we conclude that available RV data for 
GJ876 are fully consistent with the coplanar configuration (/ = 0). Nevertheless, 
we can obtain some informative upper limit on the possible non-coplanarity from 

6 To make these definitions more rigorous, we must also mention that in the Jacobi coordi- 
nate system, that we adopt here, different osculating orbits have different reference points, so 
they are no longer confocal, and this is not reflected in Fig. [8] These small displacements do 
not have practical significance, however. 
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Fig. 9 The confidence regions for the inclinational parameters (I cos<P, I sin<P), plotted in the 
same way as in Fig. See text for the detailed discussion. 



Fig. [9] the angle I likely cannot exceed 15°. However, because of the significant 
prolateness of the error ellipses in Fig. [9l the upper limit on I significantly depends 
on the value of i.e. on the orientation of the orbital nodes for planets b and c. 
So large values of I as 10° — 15° can be accepted only if the corresponding orbital 
planes intersect each other rather close to the sky tangent plane (i.e., <P close to 

or 180°). When the intersection nodes are far from the sky plane, the mutual 
inclination I is unlikely to exceed ~ 5° . 

From the non-coplanar three-planet fit bv lCorreia et al we find / = 1.9° 

and $ = 243° from their Table 2 and I = 3.5°, $ = 200° from their Table 3. 
iRivera et al. (2010) only mentioned that their non-coplanar four-planet fit yields 

1 = 3.7°. Both groups agree that there is no significant mutual inclination between 
the planets b and c. Our results do not contradict to such estimations. We have to 
admit, however, that both works discu ssed the non-coplana rity issue very briefly, 
and they omit exact uncertainties for /. ICorreia et al.l (|2010l) gave some uncertain- 
ties for ii,i2,Af2 (about 1 — 2°), but did not supply the necessary correlations, 
disabling us t o derive the inferred unc ertainties for I and/or <P. 

Previouslv. lBean fc Seifahrtl (|2009l) also tried to determine the orb it non-coplan arity 
betw een the planets b and c, based on the old Keck RV dat a from (jRivera et al.l 
2005) and HST astrometry data from (jBenedict et all 12002^ . Their non-coplanar 
orbital fit corresponds to I = 4.5° and $ = 330° in our notation. These val- 
ues also agree with co nfidence regions in Fig. [9] We must note that although 
iBean fc S eifahrt (2009) utilized the astrometic measurements, they, however, took 
into account neither the annual systematic errors in the old Keck data (Fig. [TJ), 
nor the very existence of the fourth planet, nor the red noise in RV data. They also 
do not mention whether they removed the secular acceleration effect from the RV 
data they used. Therefore, their results cannot be directl y compared with ours . 
Due to this, the possible effects from the astrometry data bv lBenedict et all ( 2002() 
remain not fully clear. We do not expect, however, that these effects are large. In- 
deed, note that the est i matio ns of, for instance, the absolute node longitudes c 
from ( Bean fc Seifahrt, 2009) possess rather large uncertainties of ~ 8° . Since the 
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astrometry data are responsible for these large uncertainties almost exclusively, 
we may suppose that such data should not constrain very much the values of 
/ ~ 5° — 10°, and probably most of the non-coplanarity information now comes 
from radial velociti es anyway. The role of the astrometric data was nevertheless 
more important in ( Bean fc SeifahrtL l2009h . since they used old RV data, which 
allowed much larger / of ~ 15° — 20°. We do not use th e astrometric data here , 
because their error properties are in fact poorly assessed. iBean fc Seifahrtl (2009) 
noted that the actual scattering of the astrometric residuals is considerably differ- 
ent from the stated instrumental uncertainties. Also, these data can easily contain 
a correlated component. We think these effects are difficult to estimate reliably, 
due to the relatively small size of the astrometric dataset. 

In the first, purely white, case reflected in Fig. O the results of the bootstrap 
simulations (not shown) are in good agreement with the confidence regions plotted, 
just like in Fig. [3] In the second, white+red, case, the bootstrap simulations are 
not helpful, since they destroy any noise correlation effects anyway. However, it is 
very likely that the formal statistical reliablity of the confidence regions in this case 
should be so high as for the white noise model. The only remaining question is how 
well the particular noise model, used in the algorithm from the Appendix lAl allows 
to eliminate the effects coming from the RV noise correlateness. To check this, we 
replotted the confidence contours from Fig. [9] assuming a bit with different noise 
models. First, we checked the picture remains practically the same for different 
simple noise correlation functions (e^' 2 ', e~ x / 2 , and 1/(1 + x 2 )). After that, we 
probed different splitting of the red part of the RV noise between the Keck and 
HARPS data. We checked the cases when the HARPS noise is purely white, and 
when it has its own red component, not tied to the Keck one. In these cases, 
the changes in the confidence regions are larger, roughly similar to the difference 
between left and right panels in Fig. [9] This may indicate that the effect of the 
red RV noise still may be not taken into account in full. It seems from the current 
data, that the HARPS red jitter may be a bit smaller than the Keck one. We 
cannot use a more accurate noise model, however. The Keck red RV jitter alone 
looks pretty estimatable without the HARPS data, but the HARPS red jitter 
is, on contrary, poorly separable. We have to live with that problem until more 
HARPS data are acquired. Anyway, all noise models tested so far, do not imply a 
significant non-coplanarity of the system and place almost the same limits on it. 



8 Long-term dynamics 

Now we can investigate the long-term dynamical evolution of the planetary system. 
We choose two best fitting orbital configurations, corresponding to the osculating 
e e = (solution I) and e e = 0.12, u) e = 45° (solution II). Both configurations are 
obtained assuming the white+red RV noise model. We tracked the evolution of 
each orbital configuration over 1 Myr term (the innermost planet d was not taken 
into account during the integration). The relative energy error was about 10 -8 
in the first integration, and about 10~ 7 in the second one. The energy error was 
walking around these values all the time and did not show a notable accumulation 
effect. This is exactly what we should expect from the symplectic integrator. We 
must note, however, that our dynamical analysis here is still rather preliminary 
and in future it is better to perform the integration over longer terms and to 
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take into account the short-period planet d using, e.g., an averaged Hamiltonian 
method (jFarago et all l2009h . 



Both orbital configurations appeared stable during the integration time. The 
evolution of the parameters related to the massive planets b and c is fairly regular 
and is close to the apsidal corotation resonance state. The main secular period of 
the system - the period of the pericenters revolution - is close to 14.3 yrs in the 
both cases. The perturbations from the fourth planet are rather small, although 
they add some minor chaotic component in the motion of the massive planets. The 
long-term evolution of this fourth planet itself looks, on contrary, considerably 
chaotic, in terms of its orbital eccentricity at least. This eccentricity, however, 
remained bounded from the upper side by 0.09 (solution I) and 0.16 (solution II). 

We also considered the evolution of the following resonant variables, corre- 
sponding to the individual two-planet resonances: 

s cb,c = 2 ^b -tpc- W C , s cb,b = 2i>b -Ipc- U b , 
s be,b = 2 ^e -1p b -U b , S bee = 2lp e - 1p b - U e , 
Sce,c = (4^ e - Vc)/3 - W c , Sce,e = (ilpe — ^c)/3 - W e . (14) 

These definitions of r esonant angles are derived from the general definitions from 
(jBeauee et all l2003h . Since we consider only coplanar configurations here, the 
angles ip and lo are replaceable by A and no (as we explained in Sect. 12. ip . The 
first pair in (|14[) corresponds to the 2/1 MMR between the planets b and c, the 
second pair - to the same resonance between b and e, and the third pair - to 
the 4/1 MMR betw een c and e. All eleven two-planet resonant angles studied by 
IRivera et al.l (|2010|) can be expressed via the six variables (1141) . Namely, 



^Pcb,c ^c6,C5 tycb : b ^cb,b: ^Pcb &cb,b ^c6,o 

fbe,b = — s be,b, <£&e,e = ~ s be,e, <Pbe = s be,b ~ s be,e, 
V-'ceO = 3Sce,c, V^ce3 = 3Sce,e, <^ce = Sce,c <Sce,e, 

^Pcel = 2Sce,c Sce,e, ty?ce2 = <Sce,c 2Sce,e- (1 5) 

We prefer to limit ourselves to the minimum possible number of variables. With 
no loss of information, we may consider only the behavior of the variables in (|14p . 

All of the quantities in (|14[) are 27r-periodic. It may seem that the quantities 
from the last pair of t| 14|) are 27r/3-periodic, but this is not strictly true. Adding 
2-7T to the angle ip c , for instance, changes s ce ,c by 27r/3 indeed, but the value of 
s C e,e is changed synchronously. It is better to say that this pair has a secondary 
vectorial period of (27r/3, 27r/3), in addition to the usual scalar period of 2ty. 

The evolution of these resonant arguments is illustrated in Fig. [10] Apparently, 
all three pairs librate in some limited regions. However, the argument s be e actually 
circulates, systematically avoiding the values |sf, ee | < 70°. The pair (s ce b , s C e,e) 
behaves in a similar manner: the argument s C e,e circulates, but usually stays in 
the range ±100° around 60° (or around 180°, or —60°, which are equivalent by 
periodicity) . For th e solution I, these results basically agree with the results by 
IRivera et al.l (l2010h . For the solution II, the libration ranges become wider, al- 



though the picture still remains qualitatively the same. 

We investigated t he beh avior of the Laplace resonance critical angle, also stud- 
ied bv IRivera et al.l (2010), sl = ipc — 3ip b + 2tp e . For the solution I, this angle 



librates between approximately ±40°, while for the solution II the relevant ampli- 
tude increases roughly twice. 
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Fig. 10 Temporal evolution (during the first 10 6 yrs) of resonant angles 114> in the GJ876 
planetary system. Each panel shows several dot-filled domains corresponding to three pairs 
of the critical angles, referring to the marked two-planet resonances. Both graphs should also 
contain one more duplicate "e:c" (green) spot, centered at (0°, 180°), but this spot is obscured 
by the "e:b" (blue) one. We removed this central e:c spot for the clarity of the figure. The third 
duplicate e:c spot (second shown) is cut half-and-half in the right part of the both graphs. See 
text for the details and discussion. 



9 Conclusions 

The orbit estimations for the GJ876 planetary system are affected by the fine effect 
of the correlated data errors, which is rather new to RV planet search surveys. The 
red RV noise was responsible, for instance, for the overestimated non-zero values 
of the planet GJ876 d eccentricity from previous works. Although we still cannot 
retract as large values of as 0.15, the available RV data are consistent with a 
circular planet d orbit. Another notable red noise effect is a systematic underesti- 
mation of the parameters uncertainties, which occure while we use the traditional 
white noise model (compare Tables [2] [3]) . For GJ876, this underestimation is typ- 
ically about 10 — 30%, and it is not removed even by rather trusted numerical 
simulation method like the bootstrap Monte Carlo. Although the correlated RV 
noise did not produce more breaking changes in the GJ876 orbital fit, the very 
existence of such type of RV measurements warns us that for other star, especially 
for those ones where the stellar jitter dominates in the total RV error budget, the 
effect of the correlated noise may be crucial. 

A simulated red noise example in Fig. [11] actually does not look like a pure 
noise at all. It looks like a bit noisy mixture of apparently periodic or maybe non- 
periodic but non-random components. If we act within the traditional white noise 
framework, these false variations may impose a huge misleading effect, and can 
ultimately lead to false planet detections. 

The planetary perturbations in GJ876 helped us to constrain not only the 
system inclination (and therefore the true planet masses), but also the planet e 
orbital eccentricity (e e < 0.2), which otherwise could not be limited better then 
by ~ 0.4. However, on contrary with the inclination, the help with the eccentricity 
constraint has a cost. The information constraining the value of e e comes indirectly 
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time 

Fig. 11 Simulated examples of the red and white noise of zero mean and unit variance. The 
white noise (background graph) looks like a thick detail-free horizontal band, whereas the red 
noise (foreground graph) shows clearly detailed structure, which leaves a very stable impression 
that such data contain some mixture of non-random periodic signals. This oscillating structure 
gets eventually suppressed at larger time intervals, but nevertheless is obvious up to tens of 
the correlation timescales (up to months in case of GJ876). The red noise was modeled as a 
Gaussian process with autocorrelation function e~ 

from irregular (probably short-term) perturbations. The irregular nature of the 
perturbations made this parameter practically indeterminable inside its admissible 
region. 

We gave improved estimations of all orbital parameters of the system, tak- 
ing into account the red noise effect and the bad determinability of the planet e 
eccentric parameters (e e ,w e ). 

We found the signs of a shallow long-term RV trend of ~ 0.2 m/(s-yr) in the 
data for GJ876. More detailed investigation suggests that the significance of the 
trend is tied to the uncertainty in the planet e eccentricity. Thus we do not claim 
that this trend is real indeed, but we believe it is an issue that should be addressed 
by future observations. If this RV trend will be confirmed, it may indicate the 
existence of an unseen distant satellite in the system. Its period should exceed 
~ 10 — 20 yrs (observational time span), and the semi-major axis should exceed, 
consequently, ~ 3 — 5 AU. At the distance from Sun of 4.7 pc, its sky separation 
should be ~ l" or more, and therefore such object could represent a good target 
for direct imaging, if it is large and bright enough. Its mass is poorly constrained, 
however: it may be as small as ~ 0.04M j up for P = 20 yrs or arbitrarily larger for 
longer P. This mass scales as oc a 2 (due to the distance-velocity law for a circular 
orbit), so for a < 100 AU that should be a planet rather than a brown dwarf. 

Finally, we investigated the system non-coplanarity, and found that although 
the current RV data are consistent with the coplanar solution, the actual mu- 
tual orbital inclination between the planets b and c should not exceed 5° — 15°, 
depending on the orientation of the corresponding mutual orbital nodes. 
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A Fitting orbits with correlated RV noise 

Strictly speaking, we have no goal here to characterize the coloured noise itself; rather, we 
want to suppress its influence on the orbital fits. In such case, we may not to care very much 
about how well our red noise models describe what really occurs. For example, it is not very 
important in practice, to which of the two "hills", the low-frequency or the unit- frequency 
one, in Fig. [2] the real correlated noise corresponds, and which of them is the alias. Both 
interpretations can almost equally explain the actual RV data we have, and therefore both types 
of the noise would eventually lead to similar effects during the RV curve fitting. Although we 
have shown in Sect. [4] that the noise non- whiteness significantly exceeds the natural "white" 
random fluctuations in the data, the latter fluctuations still largely contaminate the noise 
spectrum. It is hardly possible that small perturbations of the combined noise power spectrum 
could significantly affect our orbital fits, at least in comparison with the usual white-noise 
uncertainties. Therefore, we may limit ourselves to only very simple models of the red noise 
spectrum, which would only reflect its general monotonic decrease in frequency (which is 
actually the only robustly detectable noise spectrum property) . The secondary spectrum excess 
around the unit frequency will be reproduced automatically due to the aliasing effect in the 
time series. We also assume that the red noise component is shared between both Keck and 
HARPS data (so they differ only in the white noise parameters). It is likely that the red noise 
is just not responsible for the remaining irregular variations in the HARPS pcriodogram, since 
such variations could be expected from the white noise only. 

In addition, we may fix the shape of the autocorrelation function to some mathematically 
convenient function, which should only provide a suitably behaving frequency spectrum. In 
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this paper, we assume that the autocorrelation function of the red error noise e Ied (t) looks like 

Corr(e rcd (t), e Icd (t + At)) = p(At/r) = exp(- \At\/r), (16) 

where r is some unknown parameter characterizing the correlation timescale. According to the 
Wiener-Khinchin theorem, the corresponding noise frequency spectrum, P(f), represents the 
Fourier transform of I I16H and is equal to 

oo oo 

P(/)= I p(At/r)e 2 ^ At dAt = r J p(x)e 2 ^ dx = ^ - £ . (17) 

— oo — oo 

This function can be in good agreement with what we see in the Keck periodogram, for some 
value of the parameter r. The latter parameter characterizes the width of the low-frequency 
band that we see in the periodogram. We also considered some other realistic shapes of the 
correlation functions as alternatives: p(x) = exp(— x 2 /2), p(x) = 1/(1 + x 2 ). They produce a 
bit different frequency spectra, which are still generally similar to 1171 1. Note that the white 
noise is still present in the data and may be responsible for some constant level in the observed 
frequency spectrum. The ratio between the red and white noise contributions is a priori un- 
known. Varying r and the fraction of the red noise, we may construct a combined model power 
spectrum to be in good agreement with the observed one, even for apparently quite different 
autocorrelation functions. This means that noise parameters that we will derive from the RV 
data may be severely model-dependent and mutually correlated, but also this means that the 
red noise effect on the planetary system orbital estimations should be, on contrary, relatively 
independent on the noise correlation model (in comparison with the statistical uncertainties). 

The combined covariance function R(t±,t2) of two different RV measurements, taken at 
ti and t2 , is determined by the red noise only: 

R(t 1 ,t 2 ) = a 2 cd p((t 2 ~t 1 )/r), (18) 

whereas the variance of a single observations incorporates both white and red components: 

R(t,t) = al hitc>J +a 2 cd . (19) 

Here the index j refers to the datasct, to which the mentioned observation belongs to. That is, 
we assume that the white noise contributions are different for different datasets, whereas the 
red noise component is common. Such noise separation could take place if the red RV noise 
was actually caused by some activity effects in the stellar atmosphere. 

Since our RV time series are discrete, the covariance function J 1 8 1 1 Q P determines the fol- 
lowing N X N covariance matrix V of all available RV observations: 

V = V whitc + fT r 2 cd R rcd (r). (20) 

Here the matrix V w hit e represents the usual diagonal covariance matrix of the white part of 
the noise, and the remaining term is the common red noise component. The elements of the 
matrix R red represent pairwise correlations of the corresponding RV measurements and are 
equal to p(At/r) (with different At). When <r rod = 0, we have the usual uncorrelated white 
noise model with RV jitter values of c^ hitG .. 

With the noise model i20l . we can no longer use the objective function d 1 1) to obtain 
orbital fits, since this function does not take into account the red part of the noise. The 
correct likelihood function should now take into account the correlations between different 
RV observations. Since the likelihood function represent the joint probability density of the 
observation vector (as inferred by some fixed values of the model parameters), we must know 
the joint distribution of the whole vector of measurements, not just its covariance matrix. 
Assuming this joint distribution is multivariate Gaussian (i.e. the noise represents a Gaussian 
random process), we can replace (1101) by 

i / r Tw-i_\ 
ln£ = -- MndctV+ J-iVlnv 7 ^, (21) 



where r represents the full vector of all RV residuals {rji). We can easily see that for the white- 
only noise, when V is diagonal, the expression (121ft expands to the sum over N measurements 
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in lllOII . The divisor 7 in (1216 is again needed to reduce the bias in the variance parameters. 
The function d 2 1 1) involves the following parameters describing the RV noise structure: the 
white jitters for different datasets, c w hitej ! the red jitter, <r rcd , and the correlation timescale 
of the red jitter, r. The maximization of {ST) over the new set of free parameters (the ones 
describing the RV noise and the RV curve) yields the necessary best fitting estimations, which 
take into account the red noise effect. 



